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This  paper  uses  the  variational  approach  described  by  Walker  (2006)  for  assimilation  of  data  into  the 
nearshore  spectral  wave  model  SWAN.  The  system  uses  observed  two-dimensional  spectra  from  the  inte¬ 
rior  of  the  domain  to  correct  the  prescribed  boundary  conditions  for  the  forward  model.  The  objective 
function  that  determines  the  amount  of  correction  to  be  applied  is  derived  with  the  assumption  that 
the  differences  between  observations  and  model  predictions  are  mainly  a  result  of  specification  of  incor¬ 
rect  spectra  at  the  boundary.  Using  synthetic  data,  we  show  that  the  system  reproduces  the  correct  wave 
spectra  at  the  boundary  and  converges  to  the  solution  with  accuracy  greater  than  95%  in  only  a  few  model 
iterations.  Use  of  the  assimilation  system  to  estimate  the  wave  field  is  demonstrated  for  Santa  Rosa 
Island,  FL.  Results  show  excellent  agreement  with  independent  observations  of  the  bulk  (or  integrated) 
wave  parameters  such  as  significant  wave  heights,  peak  wave  periods  and  mean  wave  directions,  and 
good  agreement  with  observations  of  the  two-dimensional  wave  spectra.  The  accuracy  of  the  system  is 
reduced  when  there  is  relatively  little  energy  at  the  assimilation  location  or  when  the  nonlinear  processes 
due  to  wind  (such  as  active  wave  growth,  nonlinear  transfer  of  energy  between  frequencies  and  direc¬ 
tions,  and  breaking)  are  dominant  in  the  region  of  interest. 
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1.  Introduction 

Accurate  modeling  of  nearshore  processes  is  critical  to  under¬ 
standing  their  impact  on  the  coastal  region.  Dominant  processes 
in  the  nearshore  are  forced  by  the  action  of  waves.  Nearshore  cir¬ 
culation  models  as  well  as  sediment  transport  models  use  local 
wave  energy  to  model  the  currents  as  well  as  the  morphological 
changes  in  the  region.  Robust,  sophisticated  models  exist  for  pre¬ 
dicting  the  waves  as  well  as  the  resulting  hydrodynamics  in  the 
nearshore  region.  However,  the  accuracy  of  those  model  results  de¬ 
pend  both  on  the  accuracy  of  the  approximations  used  in  the  mod¬ 
el  to  represent  the  physics  and  on  the  model  inputs.  For  wave 
models  in  particular,  one  of  the  critical  model  inputs  is  the  wave 
conditions  at  the  domain  boundary.  In  the  perfect  scenario,  the 
model  utilizes  data  that  are  collected  at  or  very  near  the  bound¬ 
aries  to  determine  these  inputs.  Usually,  however,  such  data  are 
seldom  available.  This  dichotomy  between  what  is  generally  avail¬ 
able  and  what  is  needed  for  the  models  forms  an  additional  signif¬ 
icant  obstacle  to  accurate  modeling  of  the  wave  conditions  in  the 
domain. 

For  typical  nearshore  studies,  the  area  of  interest  and  therefore 
the  model  domain  is  small,  say  a  few  kilometers.  In  most  cases, 
there  is  minimal  active  generation  of  waves  by  the  wind  in  the 
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modeled  region  because  of  limited  fetch.  For  nearshore  applica¬ 
tions,  the  wave  model  most  commonly  used  in  conjunction  with 
hydrodynamic  models  is  SWAN  (Ris  et  al„  1999;  Booij  et  al., 
1999).  Modeling  suites  such  as  DELFT3D  (Lesser  et  al.,  2004)  use 
SWAN  to  model  the  wave  forcing  in  the  entire  domain,  which  is 
then  used  to  model  the  currents  as  well  as  morphological  changes 
in  the  nearshore.  Directional  resolution  is  important  in  this  case 
because  small  changes  in  direction  can  lead  to  significant  changes 
in  the  magnitude  and  direction  of  the  nearshore  wave-driven  cur¬ 
rents.  In  general,  the  data  available  in  such  regions  are  provided  by 
directional  wave  buoys  or  similar  observation  systems  that  mea¬ 
sure  the  wave  energy  spectrum.  In  this  paper,  we  look  at  using 
such  data  to  estimate  or  correct  the  boundary  conditions  such  that 
the  modeled  wave  conditions  in  the  entire  computational  domain 
are  improved. 

Although  data  assimilation  is  relatively  new  for  ocean  wave 
modeling,  there  have  been  some  efforts  to  utilize  observed  data 
in  conjunction  with  model  results  to  improve  wave  predictions. 
Most  of  these  efforts  have  concentrated  on  improving  wave  fore¬ 
casting  in  the  open  ocean,  and  they  have  been  driven  in  part  by 
the  availability  of  altimeters  and  synthetic  aperture  radars  (SARs) 
on  board  satellites  and  other  airborne  systems.  These  studies  were 
conducted  in  the  context  of  global  or  large-regional  ocean  wave 
modeling  using  the  WAM  wave  model  (Wave  Model  Development 
and  Implementation,  1988;  Komen  et  al.,  1994).  Initial  attempts  at 
correcting  the  wave  field  used  optimal  interpolation  methods  in 
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which  the  observed  significant  wave  height  (and  sometimes  mean 
wave  period)  was  used  to  nudge  the  model  result  towards  the  cor¬ 
rect  solution.  Such  methods  correct  the  initial  conditions  provided 
to  the  model  for  simulations  subsequent  to  the  time  the  observa¬ 
tions  are  collected.  These  methods  do  not  account  for  the  distribu¬ 
tion  of  wave  energy  in  the  directional  and  frequency  space.  While 
effective  in  retrieving/modeling  the  overall  energy  in  the  vicinity  of 
the  observations,  using  just  the  bulk  parameters  is  not  enough  to 
correct  the  two-dimensional  spectrum.  In  the  mid-90s,  a  scheme 
was  proposed  where  optimal  interpolation  was  used  in  conjunc¬ 
tion  with  a  reduced  two-dimensional  wave  spectra  (Hasselmann 
et  al.,  1997;  Voorrips  et  al„  1997).  The  wave  spectra  is  decomposed 
into  principal  wave  systems  where  each  wave  system  is  character¬ 
ized  by  a  few  integral  parameters.  Use  of  spectral  information, 
rather  than  wave  heights  alone,  was  found  to  result  in  better  over¬ 
all  prediction  of  wave  frequency,  direction,  and  energy  in  the  low 
frequency  regime  (Aouf  et  al.,  2006).  Although  it  includes  more 
of  the  observed  quantities,  this  method  only  corrects  the  initial 
conditions  in  the  model  domain.  Improvements  to  model  results, 
although  better  than  using  only  integrated  parameters,  were  still 
restricted  to  a  region  around  the  observation  location  (Voorrips 
et  al.,  1997). 

Relatively  little  has  been  published  where  both  the  data  and 
wave  model  physics  are  included  in  the  data  assimilation  system. 
In  this  study,  we  use  the  variational  data  assimilation  system 
developed  by  Walker  (2006)  to  correct  the  boundary  conditions 
for  SWAN  given  observations  from  the  interior  of  the  domain. 
The  dominant  processes  are  assumed  to  be  shoaling  and  refraction; 
therefore,  the  data  for  assimilation  needs  to  be  outside  of  the  surf 
zone.  We  will  show  that  this  system  can  provide  very  good  results 
in  the  interior  of  the  domain.  These  results  include  integrated 
parameters  such  as  the  significant  wave  height,  peak  wave  period, 
and  wave  direction  as  well  as  the  complete  two-dimensional 
spectrum. 


The  other  two  velocities  Ca  and  C0  are  energy  propagation  velocities 
in  the  spectral  domain  caused  by  depth  and  current  variations.  They 
are  defined  in  terms  of  the  apparent  frequency  Q  (as  seen  by  a  sta¬ 
tionary  observer),  which  includes  a  Doppler  shift  induced  by  the 
ambient  current 


Q  =  a  +  k(U  cos  9  +  Vsind). 


The  spectral  propagation  velocities  are  given  by 


dQ 

In' 


fdQ 

V9x 


sin0  + 


d  Q 

dy 


cos  6 


(6) 

(7) 


The  source  term  Stot  on  the  right-hand  side  of  (1)  is  described  in  de¬ 
tail  in  Ris  et  al.  (1999)  and  includes  the  effects  of  wind  growth  (Sin) 
and  energy  transfer  in  the  spectrum  by  nonlinear  wave-wave  inter¬ 
actions  (resonant  triad  and  quartet  interactions,  Snl).  Significant 
additional  contributors  to  the  source  term  are  various  processes 
by  which  wave  energy  is  dissipated.  These  include  white-capping 
(Sds.w).  bottom  friction  (SdSi b)  and  depth-induced  breaking  (SdSwbr). 

This  set  of  equations  can  be  solved  for  the  action  spectrum 
JV(x,s,t)  for  a  domain  subject  to  appropriate  boundary  conditions. 
For  portions  of  the  wave  spectrum  with  propagation  velocities  that 
carry  energy  into  the  domain,  the  ‘incident’  wave  spectrum  must 
be  specified  on  the  boundary.  In  the  spectral  domain,  for  most 
practical  implementations,  the  spectral  density  is  required  to  van¬ 
ish  on  the  upper  and  lower  frequency  (er)  boundaries;  this  condi¬ 
tion  is  satisfied  by  locating  the  a  boundary  far  from  the  energy- 
containing  region  of  the  spectrum.  The  boundary  conditions  in  0 
are  that  the  spectrum  is  periodic.  In  addition  to  these  boundary 
and  initial  conditions,  complete  specification  of  the  mathematical 
problem  requires  the  bathymetry  fi(x)  and  current  fields  U(x,t), 
V(x,t)  to  be  prescribed. 


3.  Assimilation  methodology 


2.  The  assimilation  system 

The  SWAN  model  (Ris  et  al.,  1999;  Booij  et  al.,  1999)  is  a  near¬ 
shore  wave-action-balance  model  which  can  predict  the  evolution 
of  the  wave  spectrum  in  coastal  regions.  The  wave-action  spectral 
balance  is  expressed  as 

f+V-(CN)=^.  (1) 

N(x,s,f)  is  the  action  spectral  density  defined  as 

N(x,s,t)  =  E(x,s,t)/(7,  (2) 

where  the  vectors  x  =  (x,y)  and  s  =  (er,0)  represent  spatial  and  spec¬ 
tral  position  respectively,  E  is  the  energy  spectral  density,  0  is  the 
wave  direction,  and  the  dispersion  relation 

<t2  =  gk  tanh  kh  (3) 

relates  the  intrinsic  radian  frequency  a  and  the  wave  number 
k  =  2n/l  with  g  representing  the  gravitational  acceleration  and  h 
the  water  depth.  In  (1),  V  =  and  C  =  (Cx,Cy,Cff,Ce) 

represent  the  wave-energy  propagation  velocities  in  physical  and 
spectral  space.  The  x-  and  y-direction  components  of  the  wave- 
propagation  velocities  are  given  by 

CX  =  U  +  Cg  cos  9,  Cy  =  V  +  Cg  sin  0,  (4) 

where  U(x,  f)  and  \7(x,  t)  are  the  x  and  y  components  of  the  ambient 
current  velocities,  specified  as  inputs  to  the  problem,  and  the  wave 
group  velocity  is 

C'-lft  ta"h“('  +h™)-  (5) 


In  this  section  we  present  the  variational  data  assimilation  ap¬ 
proach  used.  The  approach  was  originally  developed  for  the  assim¬ 
ilation  of  synthetic  aperture  radar  data  and  is  described  in  Walker 
(2006).  In  the  present  study,  a  reduced  version  of  the  algorithm  has 
been  used  for  wave-spectrum  data.  For  simplicity  and  compactness 
of  presentation,  we  adopt  a  strong-constraint  approach  using  La¬ 
grange  multipliers  (e.g.,  LeDimet  and  Talagrand,  1986),  but  the 
same  result  can  be  derived  by  starting  with  a  weak-constraint  for¬ 
mulation  of  Bennett  (1992)  and  then  taking  the  strong  constraint 
limit.  The  model  inputs  being  estimated  are  penalized  in  the  objec¬ 
tive  function  to  ensure  uniqueness  (Bennett  and  Miller,  1991).  An 
objective  function  which  is  a  positive-definite  measure  of  the  dif¬ 
ference  between  a  set  of  observations  and  the  model  predictions 
is  defined,  augmented  with  the  SWAN  model  as  a  constraint.  This 
objective  function  is  minimized  by  adjusting  the  SWAN-model 
boundary  condition. 

For  the  purposes  of  this  study,  only  stationary  conditions  are 
considered,  and  model  inputs  are  the  incident  wave  spectra  along 
the  boundaries.  Also,  this  study  concentrates  on  the  application  of 
this  modeling  system  to  limited  area  domains  in  the  nearshore  re¬ 
gion.  While  spectral  wave  models  make  simplifications  in  the  rep¬ 
resentation  of  nonlinear  energy  transfers,  wave  generation  and 
dissipation,  the  main  assumption  here  is  that,  with  correct  bound¬ 
ary  conditions  (as  far  as  the  model  is  concerned),  the  wave  spec¬ 
trum  in  the  domain  can  be  modeled  accurately.  The  second 
assumption  is  that  the  data  used  to  correct  the  boundary  condi¬ 
tions  are  observations  from  outside  the  surf  zone.  Therefore,  the 
errors  in  model  results  caused  by  inaccuracies  in  the  representa¬ 
tion  of  depth-limited  breaking  are  not  present  at  these  locations. 
Hence,  these  forcing  terms  are  omitted  in  the  calculation  of  the 
objective  function  and  consequently  the  adjoint  model.  However, 
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these  terms  are  included  in  the  forward  model.  The  third  assump¬ 
tion  is  that  the  error  in  the  specified  boundary  spectrum  is  uni¬ 
formly  distributed  along  the  boundary. 

We  wish  to  estimate  the  wave  spectrum  E(x,s)  for  a  region  V, 
that  minimizes  the  error  compared  to  a  set  of  spectrum  observa¬ 
tions.  The  estimated  spectrum  will  be  obtained  from  the  SWAN 
model  (1),  which  serves  as  a  ‘constraint’  for  the  minimization  pro¬ 
cess.  The  boundary  spectrum  Eb( s)  used  for  the  SWAN  model  is  as¬ 
sumed  spatially  uniform  and  serves  as  the  ‘control  variable’  for  the 
minimization.  Here,  we  will  first  define  an  objective  function  rep¬ 
resenting  the  constrained  minimization  problem,  develop  the  ad¬ 
joint  equations,  and  present  the  approach  used  for  the  overall 
assimilation  algorithm. 

The  objective  function  J  is  defined  as  follows.  For  a  set  of  M 
observations  at  spatial  locations  x>,  the  error  variance  between 
the  predicted  and  observed  wave  spectrum  E  can  be  expressed  as 

-i  M  p  2 

/  =  w£  /  [E&^-Erfs)  ds,  (8) 

1  i-l  •'s 

where  E,  is  an  observation  of  E  at  location  x,  and  the  integration  is 
over  the  spectral  domain  5.  As  discussed  above,  the  situation  of 
interest  is  stationary  in  time  and  one  where  the  errors  propagate 
from  the  boundary  into  the  domain;  i.e.,  the  effects  of  processes 
represented  by  the  source  terms  on  the  right-hand  side  of  (1)  can 
be  safely  ignored  in  the  adjoint  model.  As  a  result,  the  SWAN  model 
can  be  reduced  to 


of  the  integrand  is  identically  zero,  so  for  the  minimum  in  J  we  re¬ 
quire  that 


-C  •  VA  =  - 


2cr 
A A 


M 

J2(E-Ei)d(x-X  0; 

i=l 


(13) 


this  is  the  governing  equation  for  A,  the  adjoint  wave  action  spec¬ 
trum.  The  first  variation  of  J  now  consists  solely  of 

SJ  =  J  iyj  ACn  df  +  2<j)(72N6|<;‘)N/,  ds,  (14) 


and  so  we  can  identify 


dj  _  1  dj 
d Eb  cr  d  Nb 


A 

o 


C  •  n  d{  +  2</>Eb, 


(15) 


which  is  the  gradient  of  J  with  respect  to  the  boundary  spectrum. 
The  gradient  depends  on  the  boundary  spectrum  Eb  and  the  integral 
of  the  product  of  the  adjoint  spectrum  A  and  the  boundary  normal 
component  of  the  spatial  advection  velocity  C  along  the  boundary 

an. 

So,  collecting  and  organizing  our  results,  the  complete  set  of 
equations  for  the  assimilation  algorithm  consists  of:  the  cost  func¬ 
tion  J 


J  = 


i  M  p  ^  -|  2  r 

[£(xi- s)  -  Ei(s) J  ds  +  <t>  Eb(s)2  ds, 

l_'j  J  S  J  s 


(16) 


V  •  (C N)  =  0.  (9) 

If  we  require  that  our  wave  field  satisfy  (9)  and  introduce  A(x,s)  (the 
adjoint  wave  action  spectrum)  as  a  Lagrange  multiplier,  our  con¬ 
strained  minimization  problem  becomes  one  of  minimizing 


which  serves  as  a  diagnostic,  used  to  determine  convergence;  the 
SWAN  model  (1),  for  stationary  conditions 


V-(C N)=S-f- 


(17) 


J  =  jR  js  | jjj  £(E-E,)26(x-x,)  +AV  ■  (CIV)  jds  dx  +  4> 


Eids. 


(10) 


where  the  first  term  in  the  first  integral  is  recognized  as  /  and  the 
second  term  is  from  the  ‘constraint’  (9).  The  second  integral  penal¬ 
izes  the  control  variable  Eb  and  is  necessary  to  ensure  a  unique  solu¬ 
tion  (Bennett  and  Miller,  1991);  <j>>  0  sets  the  relative  weights  of 
the  boundary  spectrum  and  /  in  the  minimization  result  obtained. 

This  objective  function  J  as  defined  here  is  a  functional  of  N  (or 
equivalently  £)  and  A.  We  will  determine  the  conditions  for  a  min¬ 
imum  in  (10)  by  determining  where  the  first  variations  with  re¬ 
spect  to  N  and  A  vanish.  Taking  the  first  variation  with  respect  to 
A  recovers  the  constraint  Eq.  (9).  Prior  to  taking  the  first  variation 
with  respect  to  N,  we  recast  the  equation  in  terms  of  N,  rearrange, 
and  make  use  of  the  divergence  theorem  to  get 


J 


Ni)2, 5(x  -  Xj)  - NC • VA 


+  f  [  (VAC  •  h  d{  ds  +  (j>(T 2  j  Nl  ds, 

Js  Jon  Js 


ds  dx 


ni) 


where  the  spatial  portion  of  the  second  integral  is  over  the  domain 
boundary  81 Z,  and  n  is  a  unit  normal.  The  first  variation  with  re¬ 
spect  to  N  is 


SJ 


+l{LAt  ■nd£  +  2(jx72N(,  jdJVb  ds. 


SN  ds  dx 


(12) 


where  the  two  integrals  on  the  second  line  have  been  combined. 
Since  SN  is  arbitrary,  the  first  integral  will  vanish  only  if  the  balance 


the  adjoint  SWAN  model 
~  ~  2rr  M 

C'-VA  =  -—  J2(E-Ei)S(x-Xt),  (18) 

m  i-l 

where  C  =  -C;  and  finally,  the  gradient  of  J  with  respect  to  the 
boundary  condition  Eb 

J-=  f  ^C-h  dc  +  2<!>Eb.  (19) 

9Eb  Jl7R  o 

Eqs.  (17)— (19)  are  the  Euler-Lagrange  equations,  which  define  the 
minimum  for  J. 

It  should  be  noted  that  the  adjoint  SWAN  equation  was  devel¬ 
oped  using  a  homogeneous  form  of  the  SWAN  model,  while  in 
the  assimilation  algorithm  the  complete  SWAN  model  is  used  for 
the  forward  modeling.  If  the  observation  locations  are  restricted 
to  areas  in  the  domain  where  the  errors  primarily  arise  from  the 
propagation  of  errors  at  the  boundary,  not  because  of  dissipation 
or  generation,  then  it  is  expected  that  the  assimilation  system  will 
be  able  to  correct  the  wave  field  in  the  entire  domain.  In  addition, 
using  the  full  SWAN  model  for  the  forward  model  allows  the  full 
region  to  be  calculated  as  accurately  as  possible,  in  particular  the 
shallow  region  between  the  observation  locations  and  the  shore¬ 
line  that  includes  the  surf  zone. 

The  assimilation  algorithm  proceeds  from  an  initial  guess  for  Eb 
(which  could  be  zero)  and  calculates  an  estimate  of  the  wave  spec¬ 
trum  £(x,s)  using  (17).  The  adjoint  SWAN  model  (18)  is  then  solved 
with  homogeneous  boundary  conditions  (for  inward  propagating 
waves)  and  with  the  spectrum  error  as  the  source  term.  The  adjoint 
spectrum  at  the  boundary  A(x,s),  x  e  8TZ,  and  Eb( s)  are  then  used  in 
(19)  to  calculate  the  gradient.  The  gradient  is  used  in  a  conjugate- 
gradient  minimization  scheme  to  iteratively  determine  the  Eb( s) 
which  minimizes  J. 
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4.  Application  of  the  assimilation  algorithm 

The  assimilation  system  is  tested  off  the  coast  of  Santa  Rosa  Is¬ 
land,  FL,  where  a  study  was  undertaken  during  the  latter  part  of 
January,  2009  (Edwards  et  al.,  2009).  Wave  energy  during  this  per¬ 
iod  was  mild,  with  significant  wave  heights  (Hs)  less  than  1.5  m, 
and  the  waves  entered  the  region  mostly  from  the  southeast.  Prior 
to  the  start  of  the  experiments,  the  bathymetry  was  determined 
using  single  and  multi-beam  acoustic  surveys.  The  computational 
domain  and  the  bathymetry  contours  are  shown  in  Fig.  1.  The 
water  depth  along  the  offshore  boundary  varied  from  20  to  25  m, 
and  it  decreased  gradually  shoreward.  From  the  10  m  contour  to 
the  shoreline,  there  is  very  little  alongshore  variability  in  the  water 
depths.  Multiple  buoys  were  deployed  in  this  region  during  the 
study.  These  buoys  recorded  and  transmitted  measurements  of 
directional  wave  energy  every  hour.  The  buoy  locations  in  the  do¬ 
main  are  also  shown  in  Fig.  1.  The  corresponding  water  depth  at 
each  location  and  the  duration  of  operation  for  each  buoy  are  given 
in  Table  1.  One  of  these  buoys  (Sentry  ADCP,  henceforth  SAB)  also 
recorded  the  current  profile. 

Before  looking  at  the  results  using  the  collected  data,  we  first 
look  at  some  synthetic  data  sets  generated  by  running  the  SWAN 
model.  Using  synthetic  data  allows  us  to  evaluate  the  system  inde¬ 
pendent  of  the  accuracy  of  the  SWAN  model  and  also  allows  us  to 
verify  that  the  system  can  recreate  the  boundary  conditions.  To 
generate  the  data,  SWAN  was  run  for  the  region  shown  in  Fig.  1. 
The  incident  wave  spectrum  was  a  Pierson-Moskowitz  spectrum 
as  given  by  Donelan  et  al.  (1999), 


E(cr,d)  =  <xg2(2n)  4/  5exp 


-1.25 


^sech2[/?(0  —  Bp)], 


(20) 


which  has  the  wind  speed  (U10)  and  direction  (0P)  at  10  m  above  the 
ocean  surface  and  the  directional  spreading  (p)  as  the  defining 
parameters.  In  the  above  equation, /=  2nja  is  the  wave  frequency 
and  we  use  the  default  values  of  a  =  0.0081  and/m  =  0.13g/U]O  given 
in  Alves  and  Banner  (2003).  We  use  /?  =  2  as  the  value  for  the  direc¬ 
tional  spreading  factor.  We  show  six  different  simulations  (Table  2) 
ranging  from  little  wave  activity  to  significant  wave  energy  in  the 
domain,  either  with  winds  present  in  the  domain  (i.e.,  possibility 
of  wave  generation  in  the  domain)  or  no  winds  in  the  domain  (all 
energy  coming  from  the  boundary).  Wave  spectra  obtained  at  the 


x(m) 


Fig.  1.  Computational  domain  at  Santa  Rosa  Island,  FL,  showing  the  bathymetry  in 
the  region.  The  locations  of  the  buoys  are  also  shown:  Triaxys  buoy  1  (x),  Triaxys 
buoy  2  (+),  Sentry  ADCP  buoy  (0),  Sentry  buoy  (□). 


Table  1 

Wave  buoy  types,  the  corresponding  water  depths  at  their  deployment  locations  and 
their  duration  of  operation. 


Buoy  type  (name) 

Water  depth (m) 

Duration  of  operation 

Triaxys  #1  (TA1) 

5.3 

Jan  27,  5  pm  -  Jan  31,  7  pm 

Sentry  ADCP  (SAB) 

10.2 

Jan  27,  3  pm  -  Feb  5,  9am 

Triaxys  #2  (TA2) 

17.8 

Jan  28,  1  am  -  Jan  29,  5  am 

Sentry  (SIB) 

20.6 

Jan  27,  2  pm  -  Feb  5,  9  am 

Table  2 

Wind  speed  used  for  generating  the  boundary  spectrum  for  the  different  cases  in  the 
twin  experiments  and  the  resulting  wave-averaged  parameters.  Cases  (b),  (d)  and  (f) 
are  the  same  as  cases  (a),  (c)  and  (e),  but  with  the  corresponding  U ]0  acting  in  the 
domain  as  well.  The  mean  direction  was  1 00°  (waves  from  the  east  have  direction  0° 
and  from  the  south  have  direction  90°). 


Case 

U  io(m/s) 

Hsig(m) 

Mean  period  Tm( s) 

(a) 

5 

0.56 

2.74 

(c) 

10 

2.28 

5.76 

(e) 

12 

3.3 

6.96 

location  of  Triaxys  buoy  #1  (henceforth  TA1 )  from  the  SWAN  model 
run  was  used  as  input  to  the  assimilation  system.  In  the  three  cases 
where  winds  were  present  in  the  domain,  a  constant  value  of  U10 
was  specified  throughout  the  domain.  Furthermore,  for  these  cases, 
quadruplet  interaction  and  dissipation  due  to  white-capping  were 
enabled  in  the  forward  model. 

Fig.  2  shows  the  difference  in  the  significant  wave  height  be¬ 
tween  the  results  from  the  original  SWAN  model  run  and  the  sim¬ 
ulation  using  boundary  conditions  given  by  the  assimilation 
system  for  the  three  cases  in  Table  2,  with  and  without  including 
the  effects  of  wind  in  the  domain  for  the  forward  model.  For  the 
cases  where  the  wave  energy  enters  the  domain  only  from  the 
boundary  (cases  a,  c,  and  e),  the  difference  between  the  data  and 
the  results  is  less  than  5%  everywhere  in  the  domain.  We  see  some 
degradation  in  the  results  away  from  the  assimilation  location.  The 
errors  are  largest  for  the  case  where  the  incoming  wave  energy  is 


(a) 


0  2  4  6  8 


(b) 


0  2  4  6  8 


■ 


0 


x  (km) 


(f) 


x  (km) 


Fig.  2.  Difference  between  the  actual  and  the  modeled  significant  wave  height  ( Hsig ) 
as  a  percentage  of  the  actual  wave  height  for  all  six  cases. 
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Fig.  3.  Difference  between  the  actual  and  the  modeled  mean  wave  period  (Tm)  as  a 
percentage  of  the  actual  mean  wave  period  for  all  six  cases. 


x  (km)  x  (km) 

Fig.  4.  Difference  between  the  actual  and  the  modeled  mean  wave  direction  (0m)  as 
a  percentage  of  the  actual  mean  wave  direction  for  all  six  cases. 


Time  2009  (mm/dd) 


Fig.  6.  Comparison  of  integrated  parameters  at  the  assimilation  location  (Triaxys 
buoy  1).  Solid  lines  are  data  and  dashed  lines  are  model  results. 


Time  2009  (mm/dd) 


Fig.  7.  Comparison  of  integrated  parameters  nearshore  (Sentry  ADCP  buoy).  Solid 
lines  are  data  and  dashed  lines  are  model  results. 


0  2  4  6  8  10 

Iteration  number 

Fig.  5.  The  cost  function  as  a  function  of  the  iteration  number  for  case  (f). 


Time  2009  (mm/dd) 


Fig.  8.  Comparison  of  integrated  parameters  near  the  offshore  boundary  (Sentry 
buoy).  Solid  lines  are  data  and  dashed  lines  are  model  results. 
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very  high.  For  the  cases  where  there  is  wind  in  the  domain  in  addi¬ 
tion  to  wave  energy  entering  the  domain  via  the  boundary,  the  pic¬ 
ture  is  mixed.  Where  the  incoming  wave  energy  is  small  (case  b), 
the  errors  are  larger  at  the  boundary  (close  to  4%),  then  decrease 
to  near  zero  almost  immediately,  then  increase  to  about  3%  close 
to  the  assimilation  location.  This  change  in  error  can  be  directly 
attributed  to  the  generation  of  waves  in  the  domain.  Since  the  data 
at  the  assimilation  location  have  additional  energy  due  to  the  wave 
growth  in  the  domain,  the  system  adjusts  the  energy  entering  the 
boundary  to  compensate  for  this  input  of  energy  into  the  waves. 
Thus  the  result  of  the  assimilation  shows  a  larger  amount  of  wave 
energy  entering  the  domain  than  is  actually  measured.  For  the  case 
with  the  strongest  winds  (case  f),  white-capping  and  energy 
transfer  between  the  frequencies  are  also  factors.  Even  here,  the 
maximum  difference  between  the  data  and  results  from  the  assim¬ 
ilation  system  is  less  than  7%  in  the  domain. 

Fig.  3  shows  the  difference  in  the  mean  wave  period  between 
data  and  the  assimilation  system.  Again,  we  see  that  when  there 
is  no  wind  in  the  domain,  the  assimilation  system  reproduces  the 


original  data  very  accurately  (within  2%).  When  wind  is  present, 
errors  are  still  less  than  6%  in  the  domain  even  for  the  extreme  case 
(f).  Fig.  4  shows  the  difference  in  the  mean  wave  direction.  Here  we 
see  that  the  results  between  the  data  and  the  assimilation  system 
are  nearly  identical  everywhere.  The  maximum  difference  is 
approximately  3%.  Fig.  5  shows  the  cost  function  as  a  function  of 
the  iteration  step  for  case  (f).  Each  iteration  involves  one  run  of 
the  adjoint  model  and  one  run  of  the  forward  model  with  the  up¬ 
dated  boundary  condition.  We  see  that  most  of  the  variance  in  the 
domain  is  captured  in  the  first  iteration  itself.  After  about  seven 
iterations,  the  system  has  reached  its  convergence  limit. 

The  results  shown  next  were  obtained  by  using  the  hourly 
observations  of  the  wave  spectrum  from  TA1  to  estimate  the 
SWAN  boundary  condition  for  that  time,  assuming  stationary  con¬ 
ditions  (Figs.  6-10).  For  each  measurement  reported  by  the  buoy, 
the  initial  guess  for  boundary  spectrum  was  zero,  and  the  algo¬ 
rithm  tended  to  converge  in  10-15  iterations,  or  fewer.  Estimated 
spectra  obtained  from  SWAN  are  compared  to  spectra  from  the 
nearshore  buoy  SAB,  and  the  offshore  Sentry  buoy  (henceforth 


Fig.  9.  Observations  (left  column)  and  predictions  (right  column)  of  wave  spectral  energy  density  ( m2IHzldeg )  for  Jan  28, 3  pm  when  the  difference  between  observations  and 
predictions  of  Hs  is  the  largest.  Rows  (a),  (b),  and  (c)  are  at  buoys  TA1  (assimilation  location),  SAB  (nearshore  location)  and  SIB  (offshore  location),  respectively.  For  waves 
propagating  northward,  9m  =  90°. 
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Observed 


Predicted 


Fig.  10.  Observations  (left  column)  and  predictions  (right  column)  of  wave  spectral  energy  density  ( m2IHzldeg )  for  Jan  28, 12  pm  when  the  difference  between  observations 
and  predictions  of  Hs  is  the  largest.  Rows  (a),  (b),  and  (c)  are  at  buoys  TA1  (assimilation  location),  SAB  (nearshore  location)  and  SIB  (offshore  location),  respectively.  For  waves 
propagating  northward,  Bm  =  90°. 


SIB).  The  offshore  Triaxys  buoy  #2  (TA2)  was  operational  for  only  a 
short  time  (27  observations)  compared  to  the  other  locations  (191 
observations).  The  comparisons  at  this  location  were  similar  to 
those  at  SIB  and  are  not  shown  here. 

Fig.  6  compares  the  results  of  the  assimilation  algorithm  to  the 
observed  data  at  the  assimilation  location  (TA1 )  for  the  integrated 
parameters  Hs  (significant  wave  height),  Tp  (peak  wave  period),  and 
0m  (mean  wave  direction)  for  the  entire  observation  time  period. 
(In  this  figure  and  elsewhere,  all  times  are  UTC.)  The  maximum  dif¬ 
ference  between  observations  and  predictions  of  Hs  at  this  location 
is  0.32  m,  which  occurred  on  Jan  28  at  3  pm.  The  minimum  pre¬ 
dicted  error  was  approximately  0.01  m,  which  occurred  at  a  num¬ 
ber  of  different  observation  times.  The  mean  difference  between 
the  observations  and  the  predictions  is  0.07  m.  We  see  from  the 
figure  that  the  peak  wave  periods  (Tp)  are  also  mostly  recovered. 
The  largest  errors  for  Tp  are  between  20  s  and  30  s  during  the 
periods  where  the  observations  at  this  location  indicate  that  Tp  is 
between  30  and  40  s.  This  is  obviously  an  error  in  the  measure¬ 
ments  as  wave  periods  in  this  region  are  seldom  that  high.  Also, 
as  the  following  figures  show,  the  wave  periods  observed  at  the 


other  locations  are  0(10  s),  which  is  typical  for  this  region.  One 
explanation  for  this  error  is  that  the  energy  during  these  periods 
is  low  and  the  spectra  are  broad-banded.  Either  the  buoy  is  picking 
up  the  actual  low-frequency  motion  in  this  region,  or  the  sensor 
noise,  which  is  common  at  low-frequencies,  is  corrupting  the  sig¬ 
nal.  If  we  only  include  data  up  to  Jan  30,  12  pm  (after  which  the 
wave  energy  in  the  domain  is  consistently  low),  the  mean  differ¬ 
ence  between  the  observed  and  predicted  peak  periods  is  less  than 
1  s.  Comparison  between  the  predicted  and  observed  mean  wave 
directions  is  also  presented.  These  also  show  similar  tendencies 
to  those  of  the  wave  period.  If  all  the  data  are  included,  the  mean 
and  maximum  errors  are  significantly  larger  than  if  only  data  till 
Jan  30,  12  pm  are  included.  Table  3  summarizes  the  errors  at  all 
three  locations. 

Fig.  7  shows  the  results  at  SAB,  the  other  nearshore  observation 
location.  In  general,  the  trend  is  the  same  as  at  TA1  for  the  differ¬ 
ence  between  the  observations  and  predictions,  but  average  differ¬ 
ences  are  slightly  larger  at  this  location.  The  mean  difference  in  Hs 
is  about  0.13  m,  while  the  largest  error  of  about  0.26  m  is  slightly 
lower  than  that  observed  at  the  assimilation  location.  The 
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Table  3 

Mean  and  maximum  differences  (in  parenthesis)  between  the  estimated  and 
observed  wave  parameters. 


All  data 

All  data  till  January  30,  12  pm 

Hs(m) 

TP(s ) 

An(deg) 

Hs(m) 

TP(s) 

0m(deg) 

TA1 

0.07(0.32) 

11.69(31.2) 

20(82) 

0.07(0.32) 

0.9(2.3) 

4(9) 

SAB 

0.13(0.26) 

2.46(6.68) 

53(197) 

0.12(0.26) 

1.45(2.8) 

22(60) 

SIB 

0.25(0.44) 

2.43(6.34) 

76(212) 

0.23(0.42) 

1. 4(3.9) 

59(159) 

predicted  energy  here  is  generally  lower  than  the  observed  value. 
In  addition,  the  mean  difference  in  the  peak  period  is  slightly  larger 
than  at  TA1 ,  with  the  predictions  generally  lower  than  the  observa¬ 
tions  by  an  average  of  about  2  s.  This  average  is  skewed  upwards 
because  of  the  larger  differences  during  the  times  when  the  energy 
is  very  low.  If  the  results  at  these  times  are  not  taken  into  consid¬ 
eration  when  calculating  the  differences,  they  are  closer  to  1.5  s  on 
average.  The  wave  directions  also  show  larger  errors  at  this  loca¬ 
tion  than  at  the  assimilation  location. 

Fig.  8  shows  the  results  closer  to  the  boundary  at  SIB.  The  differ¬ 
ences  are  even  larger  here,  which  is  to  be  expected  for  a  few  rea¬ 
sons.  First,  the  wind  forcing  and  associated  nonlinear  interactions 
are  omitted,  although  even  in  small  areas  with  very  benign  sea- 
states,  one  would  expect  some  wave  generation.  Second,  the  effect 
of  dissipation  (other  than  by  bottom  friction)  is  omitted  in  the  cal¬ 
culation  of  the  objective  function.  Thus,  to  reproduce  the  energy  at 
the  assimilation  location,  smaller  waves  than  those  observed  are 
required  at  the  boundary,  which  is  far  away  from  the  assimilation 
location.  The  mean  difference  in  the  significant  wave  heights  at  SIB 
is  almost  0.20  m,  nearly  double  that  at  SAB.  The  mean  difference  in 
the  peak  period  is  similar  to  that  seen  at  SAB,  but  the  mean  direc¬ 
tional  difference  is  somewhat  larger  than  that  seen  nearshore;  it 
ranges  from  10°  to  20°. 

Fig.  9  shows  the  observed  two-dimensional  spectra  at  the  dif¬ 
ferent  buoys  (left  column)  and  those  predicted  at  the  correspond¬ 
ing  locations  (right  column)  on  Jan  28,  3  pm  when  the  discrepancy 
between  the  prediction  and  observation  of  Hs  at  the  assimilation 
location  was  the  largest.  We  notice  that  the  spectrum  at  the  assim¬ 
ilation  location  (top  row)  is  very  similar  except  that  the  predicted 
spectrum  has  slightly  smaller  directional  spread  than  the  observa¬ 
tion.  At  the  other  nearshore  location  (SAB,  second  row),  the  pre¬ 
dicted  spectrum  shows  a  smaller  directional  spread  but  broader 
frequency  spread  compared  to  the  observations.  The  predicted 
directional  and  frequency  spread  is  similar  to  that  at  the  assimila¬ 
tion  location.  The  main  reason  for  this  is  likely  the  omission  of  non¬ 
linear  interaction  in  the  objective  function.  This  narrowing  (or  lack 
thereof)  is  more  substantial  at  the  offshore  location  (third  row).  It 
is  interesting  to  note  that  the  maximum  energy  at  the  offshore 
location  is  over-predicted,  while  the  total  energy  is  under-pre¬ 
dicted,  which  leads  to  a  smaller  significant  wave  height  at  this 
location  (Fig.  8). 

Fig.  10  shows  the  comparison  for  Jan  28,  12  pm  when  the  pre¬ 
diction  and  observation  of  Hs  at  the  assimilation  location  was 
among  the  smallest.  As  expected,  here  the  comparison  of  the  spec¬ 
trum  at  the  assimilation  location  is  excellent.  But  at  the  other  near¬ 
shore  location  (second  row)  the  observed  spectrum  is  very  narrow- 
banded,  and  even  though  there  is  some  energy  in  every  directional 
bin,  the  maximum  energy  is  substantially  larger  than  predicted.  In 
contrast,  the  maximum  energy  is  predicted  to  be  higher  at  the  off¬ 
shore  location.  Again  the  observations  show  more  energy  in  every 
directional  bin.  Thus,  even  if  the  wave  heights  and  other  averaged 
wave  parameters  are  recovered  by  the  assimilation  model,  it  does 
not  necessarily  indicate  that  the  energy  distribution  is  completely 
captured.  This  shows  that  even  when  errors  propagating  from  the 
boundary  dominate,  effects  of  nonlinear  interactions  cannot  be 
completely  ignored. 


5.  Conclusions 

The  use  of  a  wave  data  assimilation  system  based  on  the  adjoint 
technique  developed  by  Walker  (2006)  is  demonstrated  here.  The 
system  corrects  the  boundary  conditions  so  as  to  improve  predic¬ 
tions  of  the  wave  energy  in  the  entire  computational  domain. 
The  adjoint  model  is  derived  using  a  strong  constraint  approach 
which  assumes  that  the  errors  between  the  model  and  data  are 
due  solely  to  the  incorrect  specification  of  the  boundary  conditions. 
It  is  assumed  that  the  problem  is  dominated  by  propagation,  and 
the  source  terms  are  consequently  neglected  in  the  adjoint  model. 
The  model  does  extremely  well  in  reproducing  wave  characteris¬ 
tics  in  the  domain.  Results  from  the  twin  experiments  show  that 
the  assimilation  system  is  able  to  reproduce  most  of  the  energy 
in  the  domain  even  when  winds  and  the  associated  nonlinearities 
are  significant.  The  integrated  wave  properties  such  as  significant 
wave  heights,  peak  wave  periods  and  mean  directions  are  well 
recovered,  even  though  the  iterative  procedure  starts  out  with  no 
energy  in  the  domain.  On  the  other  hand,  some  of  the  obvious  defi¬ 
ciencies  of  the  model  related  to  the  initial  assumptions  used  in 
deriving  the  objective  function  appear  when  the  entire  two-dimen¬ 
sional  spectrum  is  compared  to  observations.  The  transfer  of  en¬ 
ergy  between  frequencies  as  the  waves  propagate  is  not  captured 
by  the  model.  Consequently  the  spectral  shapes  from  the  model 
are  broader  than  the  measurements.  Also,  if  most  of  the  energy  is 
generated  inside  the  domain  rather  than  entering  the  domain  via 
its  boundaries  (in  which  case  the  strong  constraint  assumption  is 
invalidated),  the  model  does  not  have  the  physics  to  capture  the  ef¬ 
fects.  Even  with  those  deficiencies,  the  system  itself  is  very  robust 
and  consistent.  However,  there  are  obvious  areas  that  need 
improvement.  Future  work  will  address  the  inclusion  of  nonlinear 
effects  due  to  interactions  between  wave  and  bathymetry,  wind 
generation,  and  white  capping. 
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